/**
 * 
 */
package ua.lviv.franko.solvers.old.twoDimention;

import java.util.ArrayList;
import java.util.Iterator;

import org.jzy3d.maths.Range;

import ua.lviv.franko.integration.ConstantFunction2D;
import ua.lviv.franko.integration.Gauss2D;
import ua.lviv.franko.matrix.RowType;
import ua.lviv.franko.matrix.hash.gauss.Gauss;
import ua.lviv.franko.matrix.hash.gauss.HashMatrixGaussStyle;
import ua.lviv.franko.mlsFunctions.oneDimension.weightFunctions.CubicSpline;
import ua.lviv.franko.mlsFunctions.twoDimension.BaseFunction;
import ua.lviv.franko.mlsFunctions.twoDimension.Particle;
import ua.lviv.franko.mlsFunctions.twoDimension.Values;
import ua.lviv.franko.mlsFunctions.twoDimension.weightFunctions.IWeightFunction2D;
import ua.lviv.franko.mlsFunctions.twoDimension.weightFunctions.SquareSupport;
import ua.lviv.franko.solvers.twoDimention.EquationParameters;
import ua.lviv.franko.visualisation.EFGResult2D;

/**
 * @author marko
 * 
 */
public class OptimizedSolver {

	// general data
	protected int						quadraturePointsNumber;
	protected Gauss2D					gau;
	protected EquationParameters		parameters;
	protected HashMatrixGaussStyle		mainMatrix;
	protected double[]					b;

	// // data for background integration;
	// protected ArrayList<Rectangle> recList;
	// // protected ArrayList<Node> nodeList;
	// // protected ArrayList<NodePair> edges;
	// data for base functions
	protected ArrayList<Particle>		particles;
	protected ArrayList<BaseFunction>	baseFunctions;

	public OptimizedSolver(int quadraturePointsNumber, ArrayList<Particle> particles, EquationParameters parameters, int PolinomialDegree, IWeightFunction2D weight) {
		// initialization
		this.quadraturePointsNumber = quadraturePointsNumber;
		this.gau = new Gauss2D(quadraturePointsNumber);
		this.particles = particles;
		this.parameters = parameters;

		// creation of base functions;
		this.createBaseFunctions(PolinomialDegree, weight);
		// continuation of initialization

	}

	public EFGResult2D calculate() {
		for (int i = 0; i < b.length; i++) {
			BaseFunction bFunc = baseFunctions.get(i);
			// if
			// ((Math.abs(Math.abs(bFunc.getXi().getXi())-0.5)<0.000001)||(Math.abs(Math.abs(bFunc.getXi().getYi())-0.5)<0.000001))
			// {
			if ((Math.abs(Math.abs(bFunc.getXi().getYi()) - 0.5) < 0.0000001)) {
				// for (int j = 0; j < mainMatrix.getN(); j++) {
				// mainMatrix.setValue(i, j, 0.0);
				// }
				mainMatrix.setValue(i, i, 100000.0);
				b[i] = 0.0;
			}
		}
		System.out.println("start gauss -");
		long start = System.nanoTime();
		Gauss lu = new Gauss(mainMatrix);
		System.out.println();
		double[] res = lu.calculate(b);
		System.out.println("end gauss - " + (System.nanoTime() - start));
		return new EFGResult2D(res, baseFunctions, particles);
	}

	public void generateMatrix() {
		mainMatrix = new HashMatrixGaussStyle(RowType.ARRAY_TYPE, particles.size());
		b = new double[particles.size()];

		for (int i = 0; i < b.length; i++) {
			System.out.println(i);
			BaseFunction baseFunction = baseFunctions.get(i);
			double[][] vals = this.gaussQuadratureIntegration(Math.max(baseFunction.getXi().getXi() - baseFunction.getXi().getRho(), -0.5),
					Math.min(baseFunction.getXi().getXi() + baseFunction.getXi().getRho(), 0.5), Math.max(baseFunction.getXi().getYi() - baseFunction.getXi().getRho(), -0.5),
					Math.min(baseFunction.getXi().getYi() + baseFunction.getXi().getRho(), 0.5), baseFunction);
			int pos1 = particles.indexOf(baseFunction.getXi());
			for (int j = 0; j < vals.length; j++) {
				int pos2 = particles.indexOf(baseFunction.getNodeList().get(j));
				mainMatrix.setValue(pos1, pos2, vals[j][0]);
				mainMatrix.setValue(pos2, pos1, vals[j][1]);
			}
		}

//		MatrixLUStyle.printMatrix(mainMatrix);
//		MatrixLUStyle.printMatrix(b);

	}

	private void createBaseFunctions(int PolinomialDegree, IWeightFunction2D weight) {
		this.baseFunctions = new ArrayList<BaseFunction>();
		for (Iterator<Particle> iterator = particles.iterator(); iterator.hasNext();) {
			Particle particle = iterator.next();
			BaseFunction bf = new BaseFunction(PolinomialDegree, particle, weight);
			bf.setNodeList(this.getParticleSupport(particle, weight));
			this.baseFunctions.add(bf);
			// bf.display();
		}
	}

	private double[][] gaussQuadratureIntegration(double ax, double bx, double ay, double by, BaseFunction baseFunction) {
		// first column - baseFunction*v; second column - v*baseFunction
		double[][] result = new double[baseFunction.getNodeList().size()][2];
		// its a position of main node - node of base function
		int pos = baseFunction.getNodeList().indexOf(baseFunction.getXi());

		double f = 0.0;

		for (int j = 0; j < quadraturePointsNumber; j++) {
			double temp = 0.0;
			for (int i = 0; i < quadraturePointsNumber; i++) {
				double x = (ax + bx) / 2 + (bx - ax) / 2 * gau.getX(i);
				double y = (ay + by) / 2 + (by - ay) / 2 * gau.getX(j);
				Values[] values = baseFunction.calculateAllFunctions(x, y);
				double kxx = parameters.Kxx(x, y);
				double kyy = parameters.Kyy(x, y);
				double kxy = parameters.Kxy(x, y);
				double kyx = parameters.Kyx(x, y);
				double px = parameters.Px(x, y);
				double py = parameters.Py(x, y);
				double q = parameters.Q(x, y);
				for (int k = 0; k < values.length; k++) {
					result[k][0] += gau.getW(j)
							* gau.getW(i)
							* (values[pos].valX * values[k].valX * kxx + values[pos].valY * values[k].valY * kyy + values[pos].valY * values[k].valX * kyx + values[pos].valX
									* values[k].valY * kxy + values[pos].valX * values[k].val * px + values[pos].valY * values[k].val * py + values[pos].val * values[k].val * q);

					result[k][1] += gau.getW(j)
							* gau.getW(i)
							* (values[pos].valX * values[k].valX * kxx + values[pos].valY * values[k].valY * kyy + values[pos].valX * values[k].valY * kyx + values[pos].valY
									* values[k].valX * kxy + values[pos].val * values[k].valX * px + values[pos].val * values[k].valY * py + values[pos].val * values[k].val * q);
					// result[k][0] = gau.getW(j) * gau.getW(i) * (
					// values[pos].valX*values[k].valX+values[pos].valY*values[k].valY
					// );
					//
					// result[k][1] = gau.getW(j) * gau.getW(i) * (
					// values[pos].valX*values[k].valX+values[pos].valY*values[k].valY
					// );
				}
				temp += gau.getW(i) * values[pos].val * parameters.F(x, y);
				// temp += gau.getW(i) * coefFunc.calculate(x, y);
			}
			f += temp * gau.getW(j);
		}
		f *= (bx - ax) * (by - ay) / 4.0;
		for (int i = 0; i < result.length; i++) {
			result[i][0] *= (bx - ax) * (by - ay) / 4.0;
			result[i][1] *= (bx - ax) * (by - ay) / 4.0;
		}
		// result[pos][1] = f;
		b[particles.indexOf(baseFunction.getXi())] = f;
		return result;
	}

	private ArrayList<Particle> getParticleSupport(Particle p, IWeightFunction2D weight) {
		ArrayList<Particle> list = new ArrayList<Particle>();
		for (Iterator<Particle> iterator = particles.iterator(); iterator.hasNext();) {
			Particle particle = (Particle) iterator.next();
			if (weight.intersectedSupport(p, particle)) {
				list.add(particle);
			}
		}
		return list;
	}

	/**
	 * @param args
	 */
	public static void main(String[] args) {
		ArrayList<Particle> nodes = new ArrayList<Particle>();
		// параметри області:
		double ax = -0.5, bx = 0.5, ay = -0.5, by = 0.5;

		int n = 5, m = 5;

		double hx = (bx - ax) / (double) n, hy = (by - ay) / (double) m;

		// вузли:
		double rho = 1.05 * hx;
		for (int i = 0; i <= n; i++) {
			for (int j = 0; j <= m; j++) {
				Particle p = new Particle(ax + i * hx, ay + j * hy, rho);
				nodes.add(p);
			}
		}

		// equation parameters:
		EquationParameters params = new EquationParameters();
		params.setF(new ConstantFunction2D(10));
		params.setKxx(new ConstantFunction2D(1));
		params.setKyy(new ConstantFunction2D(1));
		params.setKxy(new ConstantFunction2D(0));
		params.setKyx(new ConstantFunction2D(0));
		params.setPx(new ConstantFunction2D(0));
		params.setPy(new ConstantFunction2D(0));
		params.setQ(new ConstantFunction2D(0));

		long start = System.nanoTime();
		System.out.println("start generation - ");
		OptimizedSolver solver = new OptimizedSolver(11, nodes, params, 1, new SquareSupport(new CubicSpline()));
		solver.generateMatrix();
		System.out.println("end generation - " + (System.nanoTime() - start));
		solver.calculate().display(new Range(-0.5, 0.5), new Range(-0.5, 0.5));

	}

}
